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I. INTRODUCTION 



We describe a method to evaluate electronic forces by Path Integral Monte Carlo (PIMC). Elec- 
tronic correlations, as well as thermal effects, are included naturally in this method. For fermions, 
a restricted approach is used to avoid the "sign" problem. The PIMC force estimator is local and 
has a finite variance. We applied this method to determine the bond length of H2 and the chemical 
reaction barrier of H+H2 — >H"2+H. At low temperature, good agreement is obtained with ground 
state calculations. We studied the proton-proton interaction in an electron gas as a simple model for 
hydrogen impurities in metals. We calculated the force between the two protons at two electronic 
densities corresponding to Na (r 3 — 3.93) and Al (r 3 — 2.07) using a supercell with 38 electrons. 
' The result is compared to previous calculations. We also studied the effect of temperature on the 

proton-proton interaction. At very high temperature, our result agrees with the Debye screening of 
electrons. As temperature decreases, the Debye theory fails both because of the strong degeneracy 
of electrons and most importantly, the formation of electronic bound states around the protons. 
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Forces are a basic quantity needed in understanding microscopic systems, e.g. they are basic inputs to molecular 
dynamics simulations and to predicting the equilibrium structures. For a system of electrons and nuclei in thermal 
equilibrium, it is a very good approximation to assume that the electrons follow the motion of the nuclei adiabatically. 
The forces exerted on those nuclei due to the fast moving electrons are known as the Born-Oppenheimer (BO) 
forces. Ehrenfest first related the force to the expectation value of the gradient of the potential which led to the 
Hellmann-Feynman theorem Q . Accurate results have been obtained using this theorem within the framework of the 
local density functional theory (LDA) ||. Unlike the energy, the force is directly related to the optimized geometry 
(i.e. when F = 0) and it allows one to probe every single nucleus in the system; one can determine the forces on each 
O ■ nucleus and optimize their individual positions concurrently. In conjunction with the total energy, the force can be 
, used to help to construct the potential energy surface. With energies known at certain grid points, a more accurate 
k> ' ^ can ^ e obtained if one knows the derivatives at those points also. 

Quantum Monte Carlo (QMC) methods are capable of treating many-body effects directly, which is essential in 
cases where the electron correlations are important. The computational demand of QMC scales as N 3 Q or less where 
N is the number of particles, while other methods, which depend on a complete representation of the many-body 
wave function or density matrix, such as configuration interaction (CI), have an exponential dependence on the size 
of the system. For QMC simulations, the total energy of the system usually has a variance which is proportional to 
the size of the system, making it difficult to distinguish the contribution of a single particle and the effect of a local 
displacement of nuclear positions. Consequently, the ability to calculate forces with QMC methods provides not only 
a good test of the accuracy of the commonly used LDA calculations but also an accurate many-body approach that 
could be applied to extended systems. 

The calculation of forces with QMC methods is a long standing important problem. There has been some progress 
|^| in the calculation of forces with both variational Monte Carlo (VMC) and diffusion Monte Carlo (DMC), however, 
they are still limited to relatively small systems, for example, LiH BH |jj and CuH j8|. 

The most straight forward approach to estimate forces would be to calculate the total energy difference between 
two sets of nuclear configurations which are close together. Then the force, for example the force component in y 
direction on nucleus z, can be approximated as: 

Fi dE _ E(X + AX iy y) — E(X) 

OXiy AXiy 



1 



where X is the position of all nuclei and Xi is for nucleus i. However, due to the statistical nature of Monte Carlo, the 
energy estimation is always associated with a statistical error <je and the error for the force is then of = y2a e / 'AXi y 
for independent sampling of X and X + AX. Hence the error diverges as AXi y — > making this simple approach 
impractical. 

The Hellmann-Feynman theorem expresses the force as the expectation value of the potential gradient with respect 
to the wave function ^S: 

F = -(V\V X V{X)\V). (2) 

This approach is well suited for Hartree-Fock or LDA type calculations where the trial wave functions are the eigen- 
functions of the Hamiltonian without statistical fluctuations. However, as illustrated in Ref. ||, the variance of this 
force estimator is infinite for Coulomb systems because of the 1 jr behavior of the potential as electrons approach a 
nucleus. Consequently, it is not possible to get a reliable estimation of the force by Monte Carlo methods using this 
estimator. 

There are other analytical derivative methods for QMC. For example, one can explicitly carry out the derivative of 
the variational or diffusion Monte Carlo energy estimator ^ : 

J^dR [) 

and express the force as a function of Vx&o, Vx-Ex, and Vj 1 ! 1 . Here ^ is the trial wave function and R is electronic 
coordinate. Note $o, which is the exact ground state wave function of the system, is unknown. One either needs 
to find a good approximation to Vx&o or use further Diffusion Monte Carlo walks to calculate it. Also, it may be 
difficult to determine Vx 1 ^- These difficulties have prevented routine calculations of forces with QMC. 

In this paper, we formulate the force as the derivative of the Born-Oppenheimer free energy with respect to the 
nuclear coordinates, and then evaluate the derivative with a Path Integral Monte Carlo (PIMC) technique. The force 
estimator is local and easy to compute. It also has a much smaller variance than that of Eq. (|l|). In the following, 
we first review the basic formulation of PIMC and then show how the force is computed. We apply this method to 
the molecules H2 and H3, and find good agreement between PIMC results at low temperature with those of accurate 
ground state calculations. To demonstrate the method works in an extended system, we study the proton-proton 
interaction in an electron gas, and compare it to LDA calculations. 



II. PATH INTEGRAL AND FORCE CALCULATION 

Path Integral Monte Carlo is a powerful computational technique which is capable of simulating boson systems (^] 
exactly and fermions [fl0| accurately. Besides total energy, many properties of the system, such as pair correlation 
function, specific heat, pressure, momentum distribution, and the boson super fluid density have been calculated. In 
this paper, we show how one can compute the electronic forces with PIMC as well. 

Consider a system of N non-relativistic particles (electrons and nuclei) interacting via Coulomb potential. The 
Hamiltonian is: 

N 

* = -E*V? + £^ (4) 



i<j 



where A.; = h 2 /2m,i. We will use atomic units throughout the paper: the unit for length is bohr and the unit for 
energy is hartree. In these units, the electron charge ej = —1 and the inverse electron mass = 1/2. 

In principle, all the properties of an equilibrium system at a finite temperature T can be determined from the 
thermal density matrix, which, for a Boltzmann system, can be expressed in the position representation as: 



p(R,R';P) = (R\e^ n \R') = \ \\(A^\ t p)-^ 2 exp 



4A,/3 



■exp[-U(R,R';P)] (5) 



where R = {r\, . . . ,rjv} is the set of all particle coordinates, (3 — 1/ksT is the inverse of the temperature. The 
expectation value of an operator O is: 

(0) = J dRdR'p(R,R';j3){R\0\R') (6) 
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where Z = exp(— (3J-) = J dRp(R, R; (3) is the partition function and T is the free energy. 

The calculation of many-body density matrix at a finite temperature is done by expanding it in terms of density 
matrices at higher temperatures: 



P (R°,R M ;I3) 



dR 1 dR 2 ■ ■ ■ dR^piR , R 1 ; r) • • ■ p{R M ~\ R M ; 



(7) 



where M is the number of time slices and r = fi/M is called the time step. 

It is much easier to obtain a good approximation to the high temperature density matrix since the system behaves 
like a classical system at high temperature. The pair product approximation jf|| has been shown to give much smaller 
errors in the density matrices compared with the primitive approximation, so a much larger time step can be used. 



In this approach, one solves the exact action 



; r) for a pair of particles and uses: 



U(R,R';r)=J2Mr- 



V ■• r ij i 



;r). 



(8) 



Errors occur only when three or more particles come close to each other. In our simulation, we used a matrix squaring 
method to numerically calculate the high temperature pair action. 



Because of their relevance, here we review a few properties of the Coulomb pair action 11 1. The classical limit of 
the action is rv(r) where v(r) is the Coulomb potential between the two particles. For large or high temperature, 
one can expand the action in powers of h or r (Wigner-Kirkwood expansion) to get on the diagonal: 



U2(nj,rij-,' 
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+ 0(r 4 ). 



(9) 



At small rij, quantum effects smooth the divergence of the Coulomb potential. As a result, the pair action and its 
coordinate derivative are finite. The cusp condition of the Bloch equation gives the slope at the origin: 



lim 



du 2 (r lJ7 r' tJ ;T) 

dm 



2(X i + X j ) 



(10) 



Now let us take into account fermi statistics. Without magnetic fields, the Hamiltonian is independent of spin 
and s z is a good quantum number. One can treat identical particles with different s z as different species and apply 
quantum statistics only to electrons of the same spin. Let V be the permutation operator of particle labels. Then for 
each spin state M: 



PF (R,R';(3) =^Y,{-lTp{VR,R';(3) 



(11) 



where pf is the fermion density matrix and N is the number of electrons in that spin state. 

Because odd permutations of fermions contribute a minus sign, a direct summation as for bosons will result in an 
exceedingly low efficiency as (3 or N increase [|l(| . The Restricted Path Integral Monte Carlo (RPIMC) [|l0| solves this 
problem by only allowing paths which do not cross the nodal surface of the fermion density matrix. The nodes are 
determined by pp(R t , R*;t) = 0. R*, the reference point, is a special point on the path. It is the value of the density 
matrix with respect to the reference point that restricts the paths. If one knows the exact 3N — 1 dimensional nodes, 
the RPIMC method is exact. In practice, the nodes are not known, one introduces a trial density matrix and uses its 
nodes instead. In this paper, we use the nodal surface of non-interacting particle systems. This has been shown to give 
accurate simulations of hydrogen plasma |i~2| and liquid 3 He Jl3| . The success of such a seemingly simple restriction 
can be understood if one takes a more realistic pair product density matrix and applies the end-point approximation, 
one finds the nodal surface is exactly that of the non-interacting system |10 . The off-diagonal corrections to the nodes 
scale as 0(t 2 ), they are important only at fairly low temperatures because the leading kinetic energy contribution is 

or 1 ). 

RPIMC gives an additional contribution to the action due to the restriction on the crossing of the nodal surface 
p0| . Locally the nodes can be approximated as hyper-planes. For small time step r, the nodal action is: 



A/-1 



u N (R°,R M ;[3) 



i=0 



1 — exp 



d l d l+1 
At 



(12) 



where d l is the distance of R l to the nearest nodes. 
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Now return to the problem of computing the quantum forces: consider a system of N n nuclei and N e electrons. We 
fix the position of nuclei (no nuclear kinetic energy) and calculate the Born-Oppenheimer free energy of the electrons 
as a function of the nuclear coordinates: 

F = F{{X x ,---,X Nn };(3). (13) 
The force exerted on the nuclei in thermal equilibrium is: 

F=-V X T=^-V X Z (14) 
pZ 

where Z is the partition function and Vx only acts on nuclear coordinates. 
Differentiating the path integral expression for the partition function: 




P=~ I ■■■JdR°---R M - 1 p(R°, R 1 ; t) ■ ■ ■ p{R M ~ 1 ,R M - r) x g g^Cg^V 



yj- dU(R\R l+1 ;r) 



(15) 



(16) 



where (. . .) denotes an average over the paths, R° = R M and R l are electronic coordinates only. The nuclear 
coordinates are not explicitly written in the action because they are independent of the time slice and fixed during 
the simulation. 

We arrive at the formula for the force estimator by expanding the action in terms of sums over pairs of charged 
particles: 

\ i=0 kl UjX I 

Note here k and I run over nuclear indices also. For small r or r', both U2{r, r';r) and du^jdr (Eq. |To| ) have finite 
values. At large r, the action approaches reiej/r. There is no divergence in this force estimator and thus the error is 
finite in contrast to estimators based on the Hellman-Feynman theorem. 

Fermion nodes cause an additional contribution to the force. The electronic nodal surface can depend on the 
positions of nuclei, and so does the distance of R 1 to the nodes: 

d* = d i ({X 1 ,---,X Nn }). (18) 

When taking the derivative of the action with respect to nuclear coordinates, the change of the nodal surface will 
contribute a force on the nuclei: 

f n = -±( ""«^;;_: "> ) (19) 

(20) 





However, this term vanishes with non-interacting nodal surfaces because the nodal positions are independent of the 
nuclear coordinates. 

We now investigate how the computer time will depend on the time step. In general, the variance of the force 
estimator, <7p, is a function of the time step r. Because of the high temperature approximations introduced in a 
PIMC simulation, the force is only exact in the limit of r — > 0. Hence the behavior of the variance at small r affects 
the overall efficiency of a calculation. To understand the dependence of a\ on r, we consider a typical, yet simple 
system: a hydrogen atom. Assuming independent samples of slices on the path, the variance of the force is: 

(21) 

where r l is the relative coordinate between the electron and the proton at time step i. The second term in the above 
equation is the square of the force (zero in this case) and is independent of r, so we only need to estimate the first 
term. 
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Both the small r (Eq. |T^) and large r (Eq. |9|) limit for the action is known. We approximate the derivative of the 
action as: 

du(r,r; T ) _ J 1 |r|(r c 

dr ~\r/\r\ 2 \r\ > r c [ZZ > 

where r c = v2At, the thermal De Broglie wavelength, is the radius inside of which quantum correction is important. 
The variance is thus roughly: 

4«i/^^)(^^y (23) 



r 2 Zj ^ ' V 8X 
(xt- 1 / 2 +0(1). (24) 

The error will diverge when r — ► 0, but only very slowly, as r -1 / 4 . 

The above estimate does not take into account how quickly independent samples can be generated or how adjacent 
time slices are correlated, so we performed an empirical study of the efficiency using our PIMC code. The efficiency 
of the PIMC force calculation is defined as: 

6r = . (25) 

^ aj.PT v ' 

It measures how quickly the variance of the force, a 2 F , decreases as a function of computer time. Here P is the number 
of Monte Carlo steps and T is the computer time per step. Fig. (m shows the efficiency as a function of r for the 
PIMC simulation of a H2 molecule. We found that £p scales as r 1 . 

The only terms that contribute in Eq. ( |l7| ) are terms involving the nucleus in question and another charged 
particle. The dominant contributions are local. The force on a nucleus mostly comes from nearby electrons, hence 
the force variance is mainly due to nearby electronic paths and independent of the total number of electrons. One can 
preferentially move electrons that are near to the nucleus and thereby obtain a better overall efficiency of the force 
calculation. We performed PIMC simulations of H impurities in an electron gas (see section IV) and found a power 
law behavior of the efficiency as a function of the number of electrons N e . As shown in Fig. (|J), £p ~ N~ 2 - 8 . 

To evaluate the path average in Eq. (p"7|), one samples the path space with multi-level Metropolis method, the 
level is chosen so that the diffusion of paths in both the coordinate and the permutation space is maximized. This 
is discussed in detail in Ref. [p|JTo| . Typically a path segment of 4 ~ 16 slices is moved at the same time. The 
permutation space is sampled by cyclic exchange of the labels of three particles followed by a path move (for RPIMC, 
a two particle permutation gives a minus sign and is not allowed) . This achieves ergodic sampling of the permutation 
space. 



III. FORCE CALCULATION FOR H 2 AND H 3 



To test the above approach, we apply it to the H2 and H3 molecules, since very accurate results for these systems 
are known. PIMC is a finite temperature method and it is well known that an isolated Coulomb system at a finite 
temperature would self-ionize. To circumvent this problem, we place the molecule in a periodic cube. When an 
electron ionizes and moves out of the simulation cell from one side, it will reenter from the other side and be captured 
again by the molecule. The properties of the ground state and low excited states of the system are unaffected if the 
cell is large enough because the wave functions corresponding to these states are very small at the cell boundary. 
As a result, these states hardly feel the existence of the periodic boundary condition and the wave functions are 
unchanged. For higher excited states and continuous states which are more spread out in space, wave functions from 
adjacent cells overlap with each other. Those states and their energy spectrum are thus changed in such a way as 
to prevent the ionization. We are only concerned with the low temperature properties of the system, corresponding 
to the lowest states of the molecule. The periodic boundary condition is well suited for this purpose. The minimum 
image convention fl4|| was used in calculating the Coulomb interaction. No long range contributions from images are 
included. 

The ground state of hydrogen molecule (H2) has two electrons with opposite spin and can be simulated as dis- 
tinguishable particles. The first excited electronic state(6 3 S+) has an energy of 0.39 above the ground state at a 
proton- proton distance of a — 1.4. At a temperature of T = 0.026 hartrees (/3 = 38.4), we reproduce the H2 ground 
state energy within statistical accuracy of 0.001 hartrees. The time step dependence of the PIMC force on r is plotted 
in Fig. (||). At r < 0.2 the PIMC result is close to the zero temperature result of F = —0.031 a.u. within error 
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bars. In the following, we use a time step of r = 0.2 with 196 time slices. We choose the simulation cell to be a cubic 
box with length L = 20.0. Calculations show that for L > 20, the boundaries do not affect the electron-proton pair 
correlation function at L/2 and a convergent result with respect to cell size is obtained. 

Fig. (^) shows the force between the two protons as a function of inter-proton distance. Throughout the paper, 
we will use the convention that the force between the two protons is positive if repulsive. Very good agreement is 
obtained with the essentially exact ground state calculations. We fit the forces to a straight line and determined the 
bond length (which corresponds to F = 0) of 1.399 ± 0.004 bohrs, while the "exact" ground state value is 1.401 Jl5| . 
The slope of the line gives the force constant of 0.356 ± 0.028, in agreement with the ground state calculation of 0.37 

a 

In Fig. (|5j) is plotted the total energy around the equilibrium position from the same PIMC run. Note here, in 
order to compare with ground state energy calculations, we corrected the time step error by extrapolating to the 
t = limit. This is a constant shift of the energy and is unimportant in calculating the force. Clearly, the total 
energy is very flat, and the dependence of the energy on distance is completely dominated by the noise. The reason is 
two-fold. Firstly, the PIMC force estimator Eq. ([l7|) has a lower variance than the finite difference estimator Eq. (Q). 
Secondly, the energy is at its minimum, while the force is not and hence changes more rapidly. The extremum of the 
energy (minimum, maximum or saddle points) are physically important, thus it is important to be able to calculate 
the force at those points. 

The system of H3 does not form a stable molecule. The interest in this system comes from the need to determine the 
barrier of the chemical reaction H+H2 — > H2+H. This is one of the simplest chemical reactions. It has been found that 
the barrier occurs at a collinear configuration of H-H-H with an equal nearest-neighbor separation a\ = ai = a = 1.757 
|jT7f . For a three electron system like H 3 , fermi statistics for spin 1/2 particles are needed. Anderson has discussed 
a cancellation scheme for this system j!?]] in the ground state which has no fixed-node error. We used the RPIMC 
with non-interacting electron nodes to study the collinear H-H-H configuration with equal bond distances a. We used 
a time step of r = 0.2 with 196 time slices. Fig. shows the force on the outer proton with respect to a. The 
position of the barrier corresponds to the point where the force is zero (it is a saddle point of the potential energy on 
the ai — a.2 plane), and we determine it to be at a = 1.799 ± 0.016 at a temperature of 0.026 hartrees. Despite the 
error due to our assumed nodal restriction which does not have a nodal force and the fact that we are simulating at 
a non-zero temperature, this result compares well with ground state calculations which give a = 1.757. 



IV. H 2 IN AN ELECTRON GAS 

Hydrogen impurities in metals have attracted much interest, both experimentally and theoretically |l8|-|22|] . In the 
simplest model of this system, one replaces the surrounding metal ions by a uniform positive charge background and 
keeps only the valence electrons. The lattice structure of the ions and effects of core electrons are neglected. One 
studies the energy or equivalently the forces between two protons as a function of their separation. The electron 
screening of the protons is important particularly at low electron densities. Such a model is a good approximation 
for simple metals like Na, Mg and Al where the conduction electrons are in delocalized plane wave states. 

The presence of protons provides particular difficulty for theoretical studies. Due to the lack of core electrons in a 
hydrogen atom, the electrons see the "bare" proton charge. Consequently, linear response theory is inaccurate p3[ . 
N0skov and his coworkers studied the H-H interaction in an electron gas, particularly, for short H-H distances (a < 2 
bohrs), with a variety of methods based on the density functional formalism and the local density approximation. 
First fl8[] , N0skov solved the equivalent Dyson equation of the Green function G of the system by projecting out the 
difference AG due to the presence of H impurities onto a finite local basis set. Later, he solved the same problem with 
an Effective Medium Theory |l9| and a self-consistent LDA calculation |2(]] . Perrot J2l],0 studied the long range part 
of the interaction and found indeed, the H-H potential oscillates at large inter proton distances, a, due to the Friedel 
oscillations of the electron screening. Perrot |22j also proposed another method to calculate the short-range H-H 
interactions. He considered the molecular binding of H2 at short a to be the dominant effect, and then determined 
the correction due to the electronic density. 

All of the above approaches give generally consistent descriptions of the proton-proton interaction, however, their 
results are quantitatively different. The calculation of the same problem with the PIMC could be a good test of the 
validity of the different theoretical approaches used above. It is also a good test of the algorithm proposed here when 
applied to an extended system. There are no other methods that can easily carry out such a calculation. 

Jellium is characterized by the dimensionless Wigner sphere radius r s , defined by r s — (3/47rn) 1 / 3 with n the 
electronic density. We carried out PIMC simulations of H2 in an electron gas of densities representing Na (r s — 3.93) 
and Al (r s — 2.07) and computed the force between the two protons. The constant background of ions gives no 
contribution to this force. Ewald sums are used to compute the long range contributions from periodic images of 
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proton and electron charges. Our supercell consisted of two protons with a fixed separation of a and N e electrons in 
a cube of size L with 

N < = 1& 3 - < 26 > 

Since some electrons could be bound to the protons and occupy a relatively smaller space, the actual surrounding 
electronic density is lower than n. However, this finite-size correction decreases rapidly as L increases as shown in 
Fig. (I?]) and discussed below. In most of the following, we used a temperature of T = 1/16 hartrees with two time 
steps of: r = 0.4 and r = 0.8. The r = limit is extrapolated by fitting the result to F(t) — F(Q) + ar 2 . A typical 
calculation with 38 electrons and r = 0.4 takes about 100 CPU hours on a SGI/CRAY Origin2000 computer. 

Fig. (0) shows the convergence of the result as the cell size increases. For the lower curve (diamond), the number 
of electrons is determined by Eq. ( p6| ) . For the upper curve (circle) , two additional electrons are added to the system 
so that the surrounding electron density is higher than n. The two results should bound the result for density n from 
above and below. As L increases, the difference between the two curves quickly reduces and at L > 10, no noticeable 
size dependence is present within error bars. 

We performed calculations with N e = 38 (L = 21 for r s = 3.93 and L = 11 for r s = 2.07) which is shown by tests to 
be large enough that the finite size effect is small with respect to the error bars. In Fig. (||) and Fig. (||), we plot the 
result against previous calculations of N0skov pa] and Pcrrot p2| . The presence of an electron gas greatly reduces 
the H-H binding with respect to a free H2 molecule. At r s = 3.93~(Fig. ||), where the electron density is low, the short 
range force is still dominated by the bonding electrons and it is almost the same as that of a free H2. As a increases, 
the strong attractive force in a H2 molecule becomes substantially weaker. The equilibrium H-H position is a w 1.5, 
larger than the H2 bond length. One explanation is that the anti-bonding state of the H-H system is also partially 
filled when it is placed in an electron gas. N0skov's LDA calculation agrees quite well with the PIMC calculation. The 
PIMC extends easily to large proton separations while N0skov's method has an increasing numerical error at large a 
due to his choice of the one-centered basis functions and is limited to a < 2. For a higher density, r s = 2.07 (Fig. |^), 
the H-H interaction is completely repulsive. N0skov's calculation, though giving the same qualitative description, is 
less repulsive than our PIMC result. This may be caused by the lack of the electron-electron correlation in the LDA. 
Perrot's calculation, which always predicts the existence of an equilibrium H-H distance even for densities as high as 
r s = 2.07, does not agree with either the LDA or PIMC results. 

We also studied the temperature effects on the electronic screening for a — 2.0. At low temperatures, electrons 
are bound to the protons and provide an attractive force which overcomes the Coulomb repulsion between the two 
protons. As the temperature increases, the electrons become more energetic and it is more and more difficult to 
confine them to the vicinity of protons. At sufficient high temperature, electrons are almost free, we recover the pure 
proton-proton Coulomb force(l/a 2 ). 

In the non-degenerate limit (nA 3 <C 1 where A = h/y/m e kBT), the electronic screening can be understood by the 
Debye theory p^j, where the screened proton potential is: 

v s (r) = -exp(-r/r ) (27) 
r 

and ro = (47re 2 /3n) -1 / 2 is the Debye length. As Fig. ( fl0| ) shows, at very high temperature (T > 2 hartrees), the result 
of Debye screening agrees with that of the PIMC calculations. As temperature decreases, the Debye picture fails both 
because of the strong degeneracy of electrons and most importantly, the formation of electronic bound states around 
the protons. 



V. OUTLOOK 



In this paper, we presented a method to calculate the electronic forces with PIMC simulations. There is no trial 
function involved in PIMC which makes it much easier to simulate a variety of physical systems at different geometries. 
The force estimator is local and thus it is both easy to calculate and scales well with the system size. Temperature is 
included naturally in the PIMC. This could be a disadvantage, for example, for the study of bond breaking at very 
low temperature. But it also allows one to study the finite temperature behavior of the system, as well as comparing 
with experiments directly. Finally and most importantly, electronic correlations are included so that one could study 
systems where the correlation is strong. 

We demonstrated the effectiveness of this method by applying it both to simple molecules and an extended electron 
gas. The approximation of the fermion nodal surface is the only uncontrolled approximation in this approach. This 
approximation could be significant when bound states form. The variational path integral technique M would allow 
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one to perform truly ground state force calculations using this method, at the expense of reintroducing a trial function. 
A natural extension of this work is to study the effects of nodes that also depend on nuclear coordinates. For example, 
one could use the nodal surface resulting from a one-particle self-consistent Hartree-Fock or LDA calculations [ ^0[ . 
One could also put in a "backflow" effect by transforming paths onto "quasi-particle" coordinates. This approach 
has been found to be very successful for liquid 3 He [^5| and 2D — 3D electron gas ^(|. It would be interesting to 
develop a Monte Carlo method similar to the Car-Parrinello approach with atoms moving on the Born-Oppcnhcimcr 
potential energy determined by QMC calculations once the force calculation can be done efficiently enough. 



ACKNOWLEDGMENTS 



This work was supported by NSF grant No. DMR-94-224-96, ONR grant No. N00014-92-J-1320 and the Department 
of Physics at the University of Illinois at Urbana-Champaign. Calculations were performed at the National Center 
for Supercomputing Applications. 



[1] P. Ehrenfest, Z. Phys. 45, 455 (1927). 

[2] R. P. Feynman, Phys. Rev. 56, 340 (1939); thesis, Massachusetts Institute of Technology, (1939, unpublished). 
[3] O. H. Nielsen and R. M. Martin, Phys. Rev. B 32 3780 (1985); Phys. Rev. B 32 3792 (1985). 

[4] D. M. Ceperley and L. Mitas, Quantum Monte Carlo Methods in Chemistry, in New Methods in Computational Quantum 

Mechanics, Advances in Chemical Physics, XCIII, ed. I. Prigogine and S. A. Rice, (1996). 
[5] B. L. Hammond, W. A. Lester, Jr. and P. J. Reynolds, Monte Carlo Methods in Ab Initio Quantum Chemistry, Chapter 

8 (1994). 

[6] J. Vrbik and S. M. Rothstein, J. Chem. Phys. 96, 2071 (1991). 

[7] C. J. Umrigar, Int. J. Quan. Chem.: Quan. Chem. Sym. 23, 217 (1989) 

[8] P. Belohorec, S. M. Rothstein and J. Vrbik, J. Chem. Phys. 98, 6401 (1993) 

[9] D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995) 
[10] D. M. Ceperley, Path Integral Monte Carlo Methods for Fermions, in Monte Carlo and Molecular Dynamics of Condensed 

Matter Systems, ed. K. Binder and G. Ciccotti, Editrice Compositori, Bologna, Italy, (1996). 
[11] E. L. Pollock, Comp. Phys. Comm. 52, 49 (1988). 

[12] C. Pierleoni, D. M. Ceperley, B. Bernu and W. R. Magro, Phys. Rev. Lett. 73, 2145 (1994). 
[13] D. M. Ceperley, Phys. Rev. Lett. 69, 331 (1992). 

[14] M. P. Allen and D.J. Tildesley, Computer Simulation of Liquids, Oxford (1987). 
[15] W. Kolos and L. Wolniewicz, J. Chem. Phys. 43, 2429 (1965). 

[16] Z.W. Sun, W. A. Lester, Jr. and B. L. Hammond, J. Chem. Phys. 97, 7585 (1992). 

[17] J. B. Anderson, Exact Quantum Chemistry by Monte Carlo Methods, in Quantum Mechanical Electronic Structure 

Calculations with Chemical Accuracy, ed. S. R. Langhoff, (1995). 
[18] J. K. N0rskov, Phys. Rev. B 20, 446 (1979). 
[19] J. K. N0rskov, J. Chem. Phys. 90, 7461 (1989). 

[20] O. B. Christensen, P. D. Ditlevsen, K. W. Jacobsen, P. Stoltze, O. H. Nielsen and J. K. N0rskov, Phys. Rev. B 40, 1993 
(1989). 

[21] F. Perrot and M. Rasolt, Phys. Rev. B 27, 3273 (1983). 
[22] F. Perrot, J. Phys.: Condens. Matter 6, 431 (1994). 

[23] Z. D. Popovic, M. J. Stott, J. P. Carbotte and G. R. Piercy, Phys. Rev. B 13, 590 (1976); P. Bhattacharyya and K. S. 

Singwi, Phys. Rev. Lett. 29, 22 (1972). 
[24] W. Ebeling, W. D. Kraeft and D. Kremp, Theory of Bound States and Ionization Equilibrium in Plasmas and Solids, 

Chapter 2.1, (1976). 
[25] R. M. Panoff and J. Carlson, Phys. Rev. Letts. 62, 1130 (1989). 

[26] Y. Kwon, D. M. Ceperley and R. M. Martin, Phys. Rev. B 48, 12037 (1993); Y. Kwon, D. M. Ceperley and R. M. Martin, 
submitted to Phys. Rev. B (1998). 



8 




FIG. 1. The efficiency £f of the PIMC force calculation for a H2 molecule as a function of time step r at an inverse 
temperature (3 = 19.2 a.u. in units of bohr 2 /hartree 2 second on a SGI/CRAY Origin2000 computer. The solid line is the power 
law fit r 1 ' 4 . 
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FIG. 2. The efficiency £f of the PIMC force calculation as a function of the total number of electrons. The system consists 
of two protons and N e electrons as discussed in Fig. (fjj). The solid line is the power law fit N~ 2 ' 8 . 
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FIG. 3. The time step dependence of the PIMC force calculation for a H2 molecule at an inverse temperature of (3 = 19.2 
a.u. The proton separation is 1.5 bohrs. The zero temperature result for the force is F = —0.031 a.u. |Jl5| 
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FIG. 4. The force between two protons in a H2 molecule as a function of inter-proton distance a. The force is plotted as 
positive when it is repulsive. Diamonds with error bar are PIMC results at /3 — 38.4 with 196 time slices. Open circles are 
ground state calculation by Kolos, et al. (is}] . The solid line is a linear fit to the PIMC result. The H2 bond length is estimated 
to be 1.399 ± .004. 
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FIG. 5. The Born-Oppenheimer energy of a Eh molecule by PIMC simulations are obtained from the same PIMC runs as in 
Fig. (Q). Time step error was corrected for by extrapolating to the r = limit. The solid line is ground state calculation by 
Kolos, et. al. Q. 
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FIG. 6. The force on the outer protons in a H-H-H collinear configuration of H3 molecule as a function of the nearest neighbor 
proton-proton distance a. Diamonds with error bar are PIMC simulations at (3 = 38.4 with 196 time slices. The solid line is a 
fit to the PIMC result. 



14 



0.10 



0.09 




as, 

CD 
O 



0.08 



0.07 



0.06 




6.0 



7.0 



8.0 



9.0 
L (bohr) 



10.0 



11.0 



12.0 



FIG. 7. The dependence of the proton-proton interaction in an electron gas on the size of the simulation cell. j3 = 16.0 
and the time step r = 0.4. The proton-proton separation is 1.5 bohrs. The electronic density corresponds to r a — 2.07. The 
number of electrons in the lower curve are determined by Eq. and have slightly lower density. The upper curve has two 
more electrons in the cell. 
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FIG. 8. The proton-proton interaction in an electron gas with r a = 3.93. The diamonds with error bars are PIMC results 
at /3 = 16. Two time step values r = 0.4 and 0.8 were used and the results extrapolated to r = 0. 38 electrons were in 
the simulation cell. The solid curve is LDA calculations by N0skov and the dashed curve is by Perrot ji^]. They were 
obtained by a fifth degree polynomial fit to the binding energy data and then taking an analytic derivative. For comparison, 
the proton-proton force in a H2 molecule [JlHj is also plotted as the long dashed curve. 
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FIG. 9. The proton-proton interaction in an electron gas with r s = 2.07. The solid curve is LDA calculations by N0skov 
and the dashed curve is by Perrot |23] . Other details are as in Fig. (H) . 
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FIG. 10. Temperature dependence of the proton-proton interaction in an electron gas with r s = 3.93 at a proton separation 
of a = 2.0. Diamonds with error bars and the solid curve are PIMC results. The dashed curve is from the Debye model. The 
arrow at F — 0.25 a.u. indicates the high temperature limit. 
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